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NASA Grant No. NAG-1-1478 


Simulation of Spiral Slot Antennas on Composite Platforms 

Project Goals 

• Develop robust finite element codes for narrow slot and other antennas on doubly curved 
platforms. 

• Investigate feeding techniques for conformal slot antennas and develop improved feed models 
in the context of FEM simulations. 

• Develop antenna miniaturization techniques using high contrast low loss materials and new 
anisotropic material 

• Deliver codes capable of analyzing a wide class of conformal antenna antennas without a need 
to use sophisticated meshing packages. 

Initial Project Plan 

Year 1 

• Develop the necessary codes for the analysis of slot antennas, including spirals 

• Fabricate and measure slot spiral and examine new feeding techniques. 

• Validate new codes 

• Develop mesh truncation schemes using new ABCs and artificial absorbers 

Year 2 

• Examine miniaturization techniques for slot spirals and other printed antennas using high 
contrast and anisotropic material 

• Measure performance of miniaturization procedures 

• Develop slot antenna analysis techniques for doubly curved platforms 

• Examine effects of curvature on antenna resonance characteristics 

Year 3 

• Integrate high frequency codes and FEM codes for pattern prediction on complex platforms 

• New antenna development with multi-frequency and multi-function characteristics 

• Further code development with interactive features 


Year 1 Progress Summary 

Our year 1 progress can be characterized with four major achievement which are crucial toward the 

development of robust, easy to use antenna analysis code on doubly conformal platforms. 

Specifically, 

• A new FEM code was developed using prismatic meshes. This code is based on a new edge- 
based distorted prism and is particularly attractive for growing meshes associated with printed 
slot and patch antennas on doubly conformal platforms. It is anticipated that this technology 
will lead to interactive, simple to use codes for a large class of antenna geometries. Moreover, 
the codes can be expanded to include modeling of the circuit characteristics. The attached 


report describes the theory and validation of the new prismatic code using reference 
calculations and measured data collected at the NASA Langley facilities. The agreement 
between the measured and calculated data is indeed impressive even for the coated patch 
configuration. 

• A scheme was developed for improved feed modeling in the context of FEM. A new approach 
based on the voltage continuity condition was devised and successfully tested in modeling coax 
cables and aperture fed antennas. An important aspect of this new feed modeling approach is 
the ability to completely separate the feed and antenna mesh regions. In this manner, different 
elements can be used in each of the regions leading to substantially more improved accuracy 
and meshing simplicity. 

• A most important development this year has been the introduction of the perfectly matched 
interface (PMI) layer for truncating finite element meshes. So far the robust boundary integral 
method has been used for truncating the finite element meshes. However, this approach is not 
suitable for antennas on non-planar platforms. The PMI layer is a lossy anisotropic absorber 
with zero reflection at its interface. The Univ. of Michigan was first to make use of this 
absorber in FEM simulations with remarkable success. We expect widespread use of this 
absorber. 

• Quite unexpectedly, during year 1 (instead of year 3 as originally planned) we were able to 
interface our antenna code FEMA_CYL (for antennas on cylindrical platforms) with a standard 
high frequency code. This interface was achieved by first generating equivalent magnetic 
currents across the antenna aperture using the FEM code. These currents were employed as the 
sources in the high frequency code. As shown in the attached figures, the calculations are in 
very good agreement with measured data. 

As planned, during year-1 we also fabricated a slot spiral antenna for operation between 1.1 Ghz 
and 2.5 Ghz. Of most importance is the stripline feed used in the design. Three models of the feed 
were investigated before arriving at a successful feed design. As shown in the preliminary 
measurements, the slot spiral performed very well, having a return loss of -30dB over the entire 
design bandwidth. 

Papers Submitted for Publication 

1 . T. Ozdemir and J.L. Volakis, ‘Triangular prisms for edge-based vector finite element 
analysis,” Submitted to IEEE Trans, on MTT 

2. J. Gong and J.L. Volakis, “An efficient and accurate model of the coax cable feeding structure 
for FEM simulations,” accepted for publication in IEEE Trans, on AP 

3. J. Gong, et. al., “Performance of an anisotropic artificial absorber for truncating finite element 
meshes” submitted to IEEE Trans, on AP. 

4. J. Gong and J.L. Volakis,” Optimal Selection of a Uniaxial Artificial Absorber Layer for 
Truncating Finite Element Meshes” submitted to IEEE Microwave and Guided Wave Letters 

Papers in Preparation 

1 . J. Gong and J.L. Volakis, “ FEM modeling of slot spirals using prismatic elements” 

2. J. Gong and J.L. Volakis, “Mixed element FEM modeling of aperture coupled stacked patch 
antennas of arbitrary shape” 

3. T. Ozdemir, J.L. Volakis and M. Gilreath, “Resonance characteristics of coated patch antennas 
on cylindrical platforms.” 

4. M. Numberger and J.L. Volakis, “A hybrid FEM/High frequency approach for antenna pattern 
calculations on complex platforms” 

In all, a total of eight papers resulted out of the Year-1 and related work. Our Year-2/3 work will 
proceed as planned. Basically, we will now concentrate on antenna design and miniaturization 
techniques for developing small scale UHF and VHF antennas. 
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Triangular prisms for edge-based vector finite 
element analysis of conformal antennas 

T. Ozdemir and J. L. Volakis 


Radiation Laboratory 

Department of Electrical Engineering and Computer Science 

University of Michigan 
Ann Arbor, Michigan 48109-2122 

June 5, 1995 


Abstract 

This paper deals with the derivation and validation of edge-based shape func- 
tions for the distorted triangular prism. Although the tetrahedron is often the 
element of choice for volume tessellation, mesh generation using tetrahedra is 
cumbersome and CPU intensive. On the other hand, the distorted triangular 
prism allows for meshes which are unstructured in two dimensions and struc- 
tured in the third dimension. This leads to substantial simplifications in the 
meshing algorithm and many printed antenna and microwave circuit geome- 
tries can be easily tessellated using such a mesh. The new edge-based shape 
functions presented in this paper are validated by eigenvalue, radiation and 
scattering, and input impedance computations involving conformal antennas 
on planar and cylindrical platforms. Also presented are results showing the 
effect of platform curvature on resonance behavior of such antennas. 
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1 Introduction 


The brick and tetrahedron are popular elements in finite element analysis of 
electromagnetic problems. The first is indeed attractive because of its sim- 
plicity in constructing volume meshes whereas the tetrahedron is a highly 
adaptable, fail-safe element. It is often the element of choice for three di- 
mensional (3D) meshing but requires sophisticated and CPU-intensive mesh- 
ing packages. The distorted prism (see Figure 1) is another volume element 
which provides a compromise between the adaptibility of the tetrahedron and 
the simplicity of the brick. Basically, the distorted prism allows for unstruc- 
tured meshing (free-meshing) on a surface and structured meshing in the 
third dimension. An approach for growing prismatic meshes is illustrated in 
Figure 2 and most volumetric regions in antenna and microwave circuit anal- 
ysis can be readily tessellated using such a mesh. As seen, once the surface 
mesh at the different surface levels is constructed, the prismatic mesh can be 
generated by simply connecting the nodes between the adjacent mesh sur- 
faces. This avoids use of CPU-intensive volumetric meshing packages and in 
many cases the user/analyst can construct the mesh without even resorting 
to a surface meshing package. Examples include printed circuits and anten- 
nas on planar platforms. Moreover, because of their triangular cross-section, 
the prisms overcome modeling difficulties associated with bricks at corners 
formed by planes or edges intersecting at small angles. 

A special case of the distorted prism is the right prism which is character- 
ized by the right angles formed between the vertical arms and the triangular 
faces [1]. The top and bottom faces of the right prism are necessarily parallel 
and equal, restricting them to a limited range of applications, namely, ge- 
ometries with planar surfaces. In contrast, the distorted triangular prism is 
almost as adaptable as the tetrahedron with the exception of cone-tips which 
are not likely to occur in printed antenna and microwave circuit configura- 
tions. 

In this paper, we introduce edge-based basis functions for the most general 
distorted prisms. These prisms have non-parallel triangular faces and each of 
their three vertical edges can be arbitrarily oriented. In the following, we first 
present the edge-based shape functions and then proceed with the derivation 
of the finite element matrix. Eigenvalue, radiation and scattering, and input 
impedance computations are used to validate the new basis functions for 
prismatic elements. Also included are data showing the effect of platform 
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Figure 1 : The distorted triangular prism shown with the directions 
of the edge vectors. 


curvature on the resonance behavior of conformal patch antennas. 

2 Vector Edge-based Basis functions 

Consider the distorted prism shown in Figure 3. The prism’s top and bot- 
tom triangular faces are not necessarily parallel to each other and the three 
vertical arms are not perpendicular to the triangular faces. On our way to 
specifying a set of shape functions for the prism, we proceed with the iden- 
tification of a triangular cross-section of the prism which can be uniquely 
defined given a point ( x,y,z ) interior to the prism. The cross-section con- 
tains this point as illustrated in Figure 3. A way to specify the nodes of such 
a triangular cross-section is to use the parametric representation 

= -f* s(r, t Tjj,), i = 1,2, 3, (1) 

where r, are the nodal position vectors of the triangular cross-section (see 
Figure 4(a)) whereas r, t and are associated with the top and bottom 
triangular faces, respectively, and 0 < s < 1. Thus, for s = 0 r, specify the 
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Figure 2: Illustration of an approach for constructing prismatic meshes (unstructured 
surface mesh and structured mesh along the third dimension) 


Perspective view Top view 



containing the interior point 

Figure 3: Variation of the triangular cross-section along the length of the prism 
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Figure 4: (a) Nodal coordinates, (b) triangular cross-section with the local 
coordinates £ and T). 


nodes of the bottom triangle and for s = 1 r, reduce to the nodes of the top 
triangle. Since s assumes values between these two limits, the cross-section 
sweeps the entire volume of the prism. Clearly, each (x,y,z) point belongs 
to one of the triangles specified by r, and consequently there is a unique 
correspondence between s and a given point in the prism’s volume. 

To determine s in terms of ( x,y,z ), we first note that (1) represents 
nine equations with ten unknowns (s and the coordinates of the three nodal 
points). Thus, additional equations are needed and these are the equations 
defining the plane coinciding with the triangle, viz. 

- + T + -+ 1 = 0 . * = 1 , 2,3 ( 2 ) 

a be 

x y z 

- + 7 + -+1 = 0 3 

a b c 

in which the constants a, b, and c correspond to the locations at which the 
plane intersects the x,y, and z axes, respectively. Substituting in (2) the 
values of y, and z, obtained from (1) yields three equations involving four 
unknowns a, b, c and s. Solving these three linear equations for i, r and i 
in terms of s and substituting them in (3), we end up with the third order 
polynomial 

Q 3 S 3 -(- <i 2^ 2 "I" OjS -J- ao = 0 (4) 

whose real root is the desired value of s. The expressions for the coefficients 
a,i can be given in closed form in terms of x,y, and z, but are rather lengthy 
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and bear no significance on the rest of the analysis. Therefore, they have 
been omitted but can easily be derived as outlined above. The appropriate 
solution of (4) in terms of a 0 ,i,2,3 is [2] 


where 



a 2 

S = (5l + S 2 ) "3^ 

(5) 

Si = 

[r + (q 3 + r 2 )^ 

(6) 

■52 = 

[r-(? 3 + r 2 # 

(7) 

r = 

7 - 3(0102 - 3a 0 a 3 ) - ^(azM ) 3 

OCI3 i 

(8) 

<1 = 

1 1 2 

~rd\ — 7G9 

3 9 2 

(9) 


and this completely specifies s in terms of x,y, and z. 

We next proceed with the derivation of the basis functions. We choose 
to represent the field variation across the triangular cross-section (defined by 
s) using the two-dimensional Whitney form [3]. A simple linear variation 
will be assumed along the length of the prism. Specifically, the vector basis 
functions for the top triangle edges can be expressed as 


Ni = d\ ( 1 / 2 ^ L 3 — L 3 V L 2 ) s 

N 2 = d 2 (L 3 VL 1 -L 1 VL 3 )s (10) 

N 3 — d 3 ( Li VZ / 2 — T 2 VTi)s 

and correspondingly those for the bottom triangle edges will be 

= d\ ( L 2 VZ, 3 — L 3 VL 2 ) (1 5) 

M 2 = d 2 (L 3 VL 1 -L l WL 3 ){l-s) (11) 

Ad3 = d 3 ( L\W L2 T 2 V L \ ) (1 s). 

The subscripts in these expressions identify the edge numbers as shown in 
Figure 1 and the distance parameters d x are equal to the side lengths of the 
triangular cross-section containing the observation point (see Figure 4(b)). 
Also, 

LiU,v) = 


6 
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Figure 5 : (a) Vector map of N 2 or M 2 (b) Vector map of Kj (c) Variation of u 
as a function of x, y, and z. 




COS Q3 


h 2 

cos a 2 
A-3 


£- 

e+ 


sin 0:3 

h 2 

sin <*2 
^3 


* 

7 * 


( 12 ) 


are the usual two dimensional scalar node-based basis functions [4] for the 
same triangle with a,- denoting the interior vertex angles and hi being the 
node heights from the opposite side. The variables £ and rj represent the 
local coordinates and are illustrated in Figure 4 (b). As required with all 
edge-based shape functions, N, • e, and M, • e, have unity amplitude on the 
ith edge whereas N, • ej = Mj • ej =0 for i ^ j. Their vector character 
is depicted in Figure 5(a) and as seen they simply "curl” around the node 
opposite to the edge on which their tangential components become unity. 

It remains to define the shape functions for the three vertical edges and 
we chose to express these by the linear representations 


Ki(£,f?) = £(£,»?) £1 (£,»?) 

K 2 (£,»?) = t)(£,7/)T 2 (13) 

K 3 (£,r?) = vU, V )LsU, V ). 

As before, L, are the node-based shape functions defined in ( 12 ) and a pic- 
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torial description of Ki is found in Figure 5(b). Of particular importance 
in (13) is the unit vector v(£,t}). It is a linear weighting of the unit vectors 
v 1 , v 2 and v 3 associated with the vertical arms (see Figure 5(c)), and is given 

by 

«(£,«?) = £ £.•(£» v)vi- ( 14 ) 

t = l 

This particular choice of v minimizes tangential field discontinuity across 
the faces. Contrary to the rectangular brick, tetrahedron and right prism, 
the edge-based vector basis functions for the distorted prism do not ensure 
tangential field continuity across the faces, nor are they divergenceless. The 
same is true for the curvilinear bricks [5] and, in the case of the distorted 
prism, the field discontinuity and divergence increase with surface curvature 
(see Figure 2). This is certainly an undesirable feature but in most cases 
(particularly when sampling at 15 or so points per linear wavelength and the 
radii of curvature are greater than a wavelength), the angular deviation of 
the vertical arms is quite small. Consequently, for all practical purposes the 
field discontinuity and divergence are negligible. 


3 Eigenvalue Computation 

In this section, we examine the validity of the presented edge-based func- 
tions. Specifically, we consider the eigenvalues of three different cavities us- 
ing the edge- based distorted prism as the tessellation element. We begin by 
first deriving the matrix elements following Galerkin’s testing. The weighted 
residuals of the vector wave equation are 


/ /X Ni# ( Vx VxE - fc o E )<^ = °> 

J j J M,- • (V x V x E — E)dV — 0, 
///v K ' • (V X V X E - hi E )dV = 0, 


i = 1,2,3, 
i = 1,2,3, 
* = 1,2,3, 


(15) 

(16) 
(17) 


in which N, , M, and K, comprise the nine edge-based vector basis functions 
defined in the previous section and E is the electric field vector. 
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The matrix equations are generated by introducing the representation 

3 

E = ]T[£, 7 vN t -(r) + £ tA ,M t (r) + £ t *K,(r)] (18) 

•=i 

where EiM and Em are the expansion coefficients, and correspond to 
the average amplitudes of the field vector at the rth edge. 

Substituting (18) into (15)-(17), and invoking the divergence theorem 
yields the element equations 

3 3 

E E jN [NNCij - klNNDij ] + E Ejm[N MQ j - k 2 Q NMD i: ) + 

j-1 j = 1 

E E jK [NKC i} - k 2 0 NKDij] = 0 (19) 

i= i 

3 3 

E Ejx[MNCij - k%MNDij] + £ E jM [MMC v - k 2 0 MMD t] ) + 
j = 1 j = l 

E E iK [MKCii - klMKDij] = 0 ( 20 ) 

j= 1 

3 3 

E E iN [KNCii - k 2 0 KND tj } + E - k 2 0 KMD tj ) + 

i= i i=i 

E EjkIKKCv - k 2 0 KKD tJ ] = 0, z = 1,2,3. (21) 


j=i 

where 

NNCu = J J J (V x N,) • (V x N/)dT (22) 

iVMC,-, = J J J y {V x NO • (V x M,)dV (23) 

WtfC',-/ = J J J [y x N;) • (V x K/)dy (24) 

MMCu = /// v (VxM t) .(Vx M,)dV (25) 

MKCu = J J J (V x Mi) • (V x K,)dV (26) 

A'A'C.y = J J J y (V xKi)»(V xK t )dV (27) 
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Region 1 ( 0 < z < z, ) Region 2 ( z,< z < z 2 ) Region 3 ( z 2 < z < z 3 ) 
Triangular cross-section Quadrilateral cross-section Triangular cross-section 





Figure 6: The three regions of integration over the volume of the distorted prism. 
Prism volume, integration subregions, subregion cross-section. 


NNDu 

= / / X N '* IWI ' 

(28) 

NMD i( 

= / / X N * • Mtdv 

(29) 

NKD it 

- 

(30) 

MMD« 

= j j J v M t »M t dV 

(31) 

MKDu 

= in M ‘' K ‘ dv 

(32) 

KKD it 

= //£*• •w* 

(33) 


The above integrals must be evaluated numerically and to integrate over the 
volume, the prism is divided into three regions (see Figure 6) of triangular or 
quadrilateral cross-section. In the case of the right prism, the integrals have 
closed form expressions and are given in the Appendix. 

Upon assembly of (19)- (21) and boundary condition enforcement, we ob- 
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(a) 


Mode 

k, cm' 1 

% Error 

(Exact) 

Prism 

Brick 

Tetra. 

TEfOi 

5.236 

-0.99 

-1.36 

0.44 

™ 110 

7.025 

-4.44 

-2.23 

0.70 

TEon 

7.531 

0.07 

-2.58 

1.00 

^201 

7.531 

-0.25 

-3.13 

-0.56 

TM„i 

8.179 

-0.31 

-2.09 

2.29 


(b) 


Figure 7: (a) Rectangular cavity discretized using right triangular prisms, 
(b) eigenvalues for the air-filled cavity. 


tain the generalized eigenvalue system 

[A]{x} = k 2 [B]{x} (34) 

in which A = k 2 represent the eigenvalues of the problem. The matrices [A] 
and [5] are square, symmetric, real and positive definite. 

Example 1: Rectangular Cavity 

The first example is the rectangular cavity shown in Figure 7(a). Results 
based on brick, tetrahedron [6] and prism discretizations are given in Fig- 
ure 7(b), where they are compared to the exact values. We observe that, 
overall, the data based on the prism are better than those based on the brick 
and, at least, as good as those associated with the tetrahedron. 

Example 2: Circular Cavity 

The tessellation of the circular cavity (drum) is shown in Figure 8(a). 
The mesh was constructed by first creating a surface grid of triangles at the 
base of the drum. The prisms were then generated by growing the mesh 
along the axis of the cylinder. In this case, the resulting volume element is 
the right prism and the eigenvalues are given in Figure 8(b) along with the 
exact. Except for the higher order mode TE? n, the remaining eigenvalues 
were computed to within three percent of the exact. 

Example 3: Pie-shell Cavity 

The third example is a pie-shell sector as shown in Figure 9(a). It is ob- 
tained by bending the rectangular cavity considered earlier and the resulting 
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Radius = 1cm 



(a) (b) 


Figure 8: (a) Circular cavity discretized using right triangular prisms, 
(b) eigenvalues for the air-filled cavity 



(a) (b) 

Figure 9: (a) Pie-shell cavity discretized using distorted triangular prisms, 
(b) eigenvalues for the air-filled cavity 
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no 



Figure 10: The geometrical set up for finite element-ABC formulation of radiation and 
scattering from cavities in infinite ground planes. 

volume element is the distorted prism with a vertical arm angular deviation 
of five degrees. The computed and exact eigenvalues for the first five domi- 
nant modes are given in Figure 9(b), and these testify to the accuracy of the 
distorted prisms in modeling curved geometries. 

4 Radiation and scattering from cavity backed 
structures 

As the next set of tests, we have looked at radiation and scattering from 
cavity backed radiators recessed in an infinite planar or cylindrical metallic 
platform. The general set up geometry is shown in Figure 10. As it is seen, a 
three dimensional cavity is discretized using the triangular prisms. The finite 
element mesh is extended outside the cavity in all directions about a fraction 
of the wave length and truncated using the socond order ABC reported in 

[7]. 

Example 1: Radiation from dipoles in a cavity on planar platform 

We have first looked at the radiation from a pair of current elements 
residing inside the cavity; one horizontally and one vertically oriented as 
shown in Figure 11. Radiation pattern is shown in Figures 12a and b. The 
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(Side view) 


y-directed element Ac,ual g riddin 8 z-directed element 


0.525 K 



0.75 K 


t y (Top view) 

I x 


Figure 1 1 : Locations of current elements with respect to the gridding inside the cavity 


solid and dotted line results were obtained using brick elements and closing 
the finite element domain at the aperture of the cavity via boundary integral 
(BI) method, which is an exact formulation and is explained in [11]. The 
triangular prism results were obtained in conjunction with the absorbing 
boundary condition (ABC) to truncate the mesh 0.3A o away from the cavity 
aperture in all directions. As the plots clearly demonstrate, the triangular 
prisms together with the ABC perform well. 

Example 2: Plane wave scattering from a cavity on planar platform 

Next, we have looked at the plane wave scattering from a larger cav- 
ity whose dimensions in x,y, and z directions are 1.2A e , 0.75A o , and 0.3A o , 
respectively. The geometrical set up for the finite element-ABC formula- 
tion is the same as in the previous radiation problem which is illustrated in 
Figure 10. Figures 13 and 14 show bistatic and backseat tering radar cross 
sections (RCS), respectively, of the cavity. It is seen that the FEM-ABC 
method using triangular prism elements exhibit sufficient accuracy in com- 
puting plane wave scattering from cavities in a ground plane. 

Example 3: Radiation from a patch antenna on cylindrical platform 

As the last test, we have looked at the radiation from a cavity backed 
patch antenna on a circular-cylindrical platform shown in Figure 15(a). The 
distorted triangular prisms used to discretize the cavity in Figure 9 have been 
used here also to discretize the computational domain. Finite element mesh 
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• A <«> <r/X J (dB) 


(a) 


(b) 



Observation angle 0 (deg) Observation angle + (deg) 


Figure 12: Radiation pattern of two current elements located inside the cavity; 
(a) <t>=45°; (b) 0=45° 



Observation angle e (deg) Observation angle ♦ (deg) 


Figure 13: Bistatic RCS of a cavity of dimensions 1.2 X 0 x 0.75 x 0.3 for 
a plane wave with its electric field vector making an angle of 45 witlPthe unit 
vector in 0 direction incident on the cavity aperutre at 0 mc =45 ° <j) ,,ic =30 Q , 

(a) <t>=60 a , (b) 0=45°. 
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o/r (dB) 



Observation angle G (deg) Observation angle + (deg) 

(a) (b) 

Figure 14: Backscatter RCS of a cavity of dimensions 1.2A, 0 x 0.75 X 0 x 0.3X o 
for an incident plane wave with its electric field vector making an angle of 30° 
with the unit vector in 0 direction; (a) <j)=30°; (b) 0=60° . 

has been terminated 0.3 wavelengths away from the cavity aperture by the 
first order cylindrical ABC derived in [7]. Figure 15(b) shows the radiation 
pattern computed as compared to the more rigorous finite element-boundary 
integral solutions [8]. Again the results are excellent. The reference solution 
(FE-BI) has been proven to be accurate by comparing against experimental 
data as shown in Figure 16(a) which displays a rectangular patch antenna 
backed by a continous wraparound cavity recessed in a circular ogive cylin- 
der [9]. Figure 16(b) shows the same antenna but this time with a dielectric 
superstrate covering the entire surface of the test body. Here the experimen- 
tal data is compared against finite element-ABC solution (FE-ABC) which 
terminates the finite element mesh 0.3 wavelengths away from the surface of 
the cylinder [10]. Agreement is again good. 

5 Input impedance computations 

As the last set of tests, we computed input impedances of a rectangular and a 
circular patch antenna backed by a cavity and recessed in an infinite metallic 
platform. The results are compared to more accurate boundary integral 
solutions for planar platforms, and the effect of curvature on the resonant 
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Figure 15: Radiation from a rectangular patch on circular-cylindrical platform: (a) geometiy, 
(b) radiation pattern. 
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Figure 16: Comparison with experimental results, (a) rectangular patch on a 
wraparound cavity, (b) with a dielectric superstrate, (c) test body: circular 
ogive cylinder. 
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frequency is observed for cylindrical platforms. 

Example 1: Rectangular patch on planar platform 

As the first test, we have looked at the input impedance of a 2cm x 
3cm rectangular patch on the aperture of a 5cm x 6cm x 0.079cm cavity 
recessed in a planar metallic ground plane (see Figure 17). The meshing 
of the computation volume is the same as in Figure 10 but this time the 
mesh is terminated using a metal- backed isotropic absorber instead of an 
absorbing boundary condition. This new termination technique has already 
been shown to work well in these types of geometries [12]. In Figure 17(a), the 
results refered to as FEMA-PRISM have been obtained using prisms as the 
discretization volume element in conjuction with the absorber to terminate 
the mesh. The results refered to as FEMA-BRICK, which are the accurate 
ones, have been obtained using bricks as the discretization volume elements 
in conjuction with the boundary integral technique to terminate the mesh 
at the cavity aperture [11]. The error in estimating the resonant frequency 
(the frequency at which the resistance is maximum and the reactance is 
zero) is about 20 MHz. This is less than 1% of the resonant frequency (3.17 
GHz), which is well within the range of accuracy of the moment method 
techniques. It does not make sense to compare the input resistance values 
at resonance since both techniques use an overly simplistic modeling of the 
feeding structure and the level of input resistance is very much dependant 
on the accuracy of the model. 

For each frequency computation, the metal wall at the back of the ab- 
sorber has been kept at 3cm (0.3A o at resonance) away from the edges of and 
above the cavity aperture as shown in Figure 17(b). 

The patch is fed 0.375cm off the center along the longer dimension so 
the shorter edges are the radiating edges as shown in Figure 17(c) and the 
radiated electric field is polarized in the direction of the longer edge. 

Example 2: Circular patch on planar platform 

As the second test, we have looked at the input impedance of a circular 
patch of radius 1.3cm on the aperture of a circular cavity of radius 3.47cm and 
depth 0.45cm recessed in a planar metallic ground plane (see Figure 18). In 
principle, the meshing of the computation volume is the same as in Figure 10 
but the surface mesh is as in Figure 8(a). As in the preceeding example, the 
mesh is terminated using a metal-backed isotropic absorber. In Figure 18(a), 


19 



or X, (Ohms) 


Input Impedance vs. Frequency 


1 1cm 



Finite Element-Boundary Integral 

(FEMA-BRICK) 

Finite Element-Metal Backed Absorber 

(FEMA-PRISM) 

Resonant frequency = 3.17 GHz 

Error in resonant frequency estimation < 1% 


. 5cm | 



£,=2.17 
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Figure 17: Rectangular patch on planar platform, (a) input impedance computations, (b) antenna 
geometry and absorber termination boundary, (c) radiating sides of the patch. 
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Figure 18: Circular Patch on Planar Platform, (a) Input impedance computations, (b) Antenna 
geometry and the absorber termination boundary, (c) Radiating sides of the patch 
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the results refered to as FEMA-TETRA, which are the more accurate ones, 
have been obtained using tetrahedra as the discretization volume elements 
in conjuction with the boundary integral technique to terminate the mesh at 
the cavity aperture [13]. The error in estimating the resonant frequency is 
about 50 MHz. This is 1% of the resonant frequency (3.51 GHz), which is 
again within the range of accuracy of the moment method techniques. 

For each frequency computation, the metal wall at the back of the ab- 
sorber has been kept at 2.6cm (0.3A o at resonance) away from the edges of 
and above the cavity aperture as shown in Figure 18(b). 

The patch is fed 0.8cm off the center of the patch so, as shown in Fig- 
ure 18(c), the glowing parts of the circular patch are two apposing sides 
symmetrically located with respect to the center of the patch at the same 
angular location as the feed. 

In case of both examples, the agreement between the absorber and bound- 
ary integral termination results improved as the distance of the absorber 
termination boundary from the cavity apereture was increased, as expected. 

Example 3: Rectangular and circular patches on circular- cylindrical plat- 
form, 

As the last example, we looked at the effect of the curvature on the reso- 
nance behavior of the patch antennas using prisms. We took the rectangular 
and circular patch antennas examined in Example 1 and Example 2, respec- 
tively, and curved the platforms circularly first along x-axis with R x = 5cm, 
Ry = oo then along y-axis with R r = oo, Ry = 5 cm where R x and Ry are the 
radii of curvature along x and y-axis, respectively. While curving the plat- 
form, linear dimensions of the cavity aperture have been kept the same as in 
the planar case. Figures 19(a) and (b) show the computed input resistance 
as a function of frequency for rectangular and circular patches, respectively. 
Each figure contains three input resistance plots showing the shift in the res- 
onant frequency as a result of curving the patches. First observation is that 
the resonant frequency increases with the curvature. Secondly, the amount 
of increase depends on the way the platform is curved with respect to the 
radiating sides of the patch: When the platform is curved in a way to bring 
the radiating sides closer to each other ( R x = oo, Ry = 5cm), the shift in 
the frequency is greater than (almost twice as much as) the case in which 
the radiating sides stay the same distance apart as the platform is curved 
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(Af) x = 70 MHz 
(Af) y = 150 MHz 


Figure 19: Shift in the resonant frequency due to curving, (a) rectangular patch, 
(b) circular patch 


23 



(R x = 5 cm, Ry = oo). This observation is true for both the rectangular and 
circular patches. 


6 Conclusions 

The distorted prism is an indispensable tool for the analysis of many dou- 
bly curved antenna and microwave geometries. It provides simple volume 
meshing without compromising accuracy (see the comparison with the tetra- 
hedron in Figure 7(b)). In this paper, we introduced a set of edge-based 
shape functions for the distorted prism and used them to generate the ele- 
ment equations. For validation purposes, the eigenvalues of three different 
cavities (rectangular, cylindrical and pie-shell) were computed and compared 
to the exact solutions. In almost all cases, the error was less than three per- 
cent. Radiation, scattering and input impedance computation results are 
also excellent when compared to the finite element-boundary integral solu- 
tions. Having validated the new prismatic elements and the basis functions, 
we used them to look at the effect of curvature on the resonance berhav- 
ior of conformal patch antennas. The result is that the resonance frequency 
increases with the curvature irrespective of the type of the patch. 


Appendix 

For the right prism, closed form expressions of the integrals in (22)- (33) are 
possible. Referring to Figure Al, we have 
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Figure Al: Right triangular prism: (a) Indexing of nodes and edges, (b) parameters 
defining the triangular cross-section. 
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I. Introduction 

It has been reported that a hybrid finite element-boundary integral technique 
can be employed for characterizing conformal antennas of arbitrary shape. 
Indeed, planar/non-planar, rectangular /non-rectangular designs, ring slot 
or spiral slot antennas with probe feed, coaxial cable feed, or microstrip 
line feed, can be simulated with this formulation. This is because of the 
geometrical adaptability of tetrahedral elements used in the implementation. 
However, in practice, certain configurations require extremely high sampling 
rate to accommodate the computational domain. Among them is a multi- 
layered geometry with thin substrate(s) or bonding film(s) in the presence of 
a thick spacer is a typical design, which is found attractive for increasing the 
bandwidth of the microstrip antenna. Another example is the thin (spiral) 
slot antenna, where the slot width is much smaller than other dimensions 
(diameter or inter— distance between spiral turns). In both cases, the mesh is 
extremely dense (with over 50 or even 100 samples per wavelength depending 
on geometry) whereas typical discritizations involve only 10-20 elements per 
wavelength. This redundant sampling problem is especially severe for 3-D 
tetrahedral meshes, where the geometrical details usually distort the tetras. 
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The numerical system assembled from this type of mesh often leads to large 
system conditions due to a degraded mesh quality. Also, the mesh generation 
is difficult and the solution CPU time is large. 

In this report, we propose a finite element-boundary integral formula- 
tion using edge-based triangular prism elements. It can be shown that this 
element choice is ideally suited for planar antenna configurations, include 
spirals, circular and triangular slots. Among the many advantages of these 
elements, the most important is the simplicity of the mesh generation. Also, 
it requires much smaller number of unknowns and efforts for accurate model- 
ing of complex geometries without significant increasing the system size. As 
is typically the case with finite element formulations, layered inhomogeneous 
or anisotropic material are readily modeled. This is particularly important 
since the purpose of our study is to investigate antenna miniaturization tech- 
niques and resonance control. 


II. Edge— Based Triangular Prism Elements 

Edge-based elements for vector field modeling are free of spurious modes 
and inherently satisfy the tangential continuity boundary conditions. In this 
section, we describe two types of shape functions for triangular prisms and 
briefly derive the corresponding FEM matrix elements. 


Shape Functions: Type One 

Consider a triangle shown in figure 1(a), where the three edges are numbered 
in the counterclockwise direction. In this numbering fashion, the relations 
of the edges and the corresponding two nodes of each edge are uniquely 
determined and tabulated as follows 


edges 

nodel 

node2 

i 


*2 

i 

© 

© 

2 

© 

© 

3 

© 

© 


Note that the surface area may be written as S e = — h, where /, is the length 

2 

of the i th edge and h is the height from this edge to the corresponding vertex. 
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Figure 1: (a) Triangle edge and node numbering system; (b) Geometrical 
interpretation of eq(2) and eq(3) for constructing shape functions. 


Z 



8 


i=1,2,3 j=4,5,6 k=7,8,9 


Figure 2: Illustration of a typical triangular prism element and the numbering 
system. 
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It is also noted that the height may be expressed as 


h = ti • z x (r - p,-), ( 1 ) 

where r is a position vector right on the i th edge, r, is the position vector at 
the vertex and e, is the unit vector along the edge. Hence, the triangle area 
becomes 


56 = \ hi ' £ x ( r “ r 0- (2) 

Referring to figure 1 (b), it is evident that if r is measured within the triangle, 
then (2) gives the area of another triangle (12'3') with the same length of the 
base edge shift upward so that r is located on it. In this case, if one writes 
a vector S, such that 


s *' = 2 ^ x ( r “ r «)> (3) 

then e • S, becomes the ratio of this triangle area to the area S e of triangle 
(123) and this is the design of the edge basis functions for (2-D) triangular 
elements. The corresponding volumetric basis functions associated with the 
triangular prism shown in figure 2 may then be written, by inspection, as 


V, 

Vi 

V* 


_ (s-A z) 

Az b ’ 

_ (zg± Az-z) 

Az 

= Kk 


i = 1,2,3 

i = 4,5,6 
lb = 7,8,9 


(4) 


where Ck is the triangle simplex coordinate associated with the vertex ( x k , y k ) 
to which the k th vertical edge is linked; z c and Az are the z coordinate of 
the base surface and the height of the prism, respectively. It should be 
remarked that the basis functions of this format are divergence- less. In order 
to simplify our formulation, we will rewrite the above basis functions in 
Cartesian coordinates, yielding, 

v ‘ = 2 _ Hy - y«)* + (* - *.-)y] *' = 1,2,3 

Vj = 2S^Az ^ Zc + Az ~ z ^ [“(2/ - yj)z + (* - Xj)y] j = 4, 5 , 6 ( 5 ) 

V* = z[(xkiVk 2 - Xk 2 Vki) - x{yk 2 - Vki) + y(x k 2 - Xki)] k = 7, 8, 9 
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where (x t , y t ) and (xj, y : ) are the vertices of the i th and j ih edges, respectively; 
(x k ,yk), {xki,Vki) and ( x k 2 ,Vk 2 ) are three nodes or vertices in counterclock- 
wise rotation. It is straightforward to derive the curl of (5) and they are 


given by 



V x V,- 

= - 2ge Az ((* *<)* + (y - yi)y - z ( z ~ 2 ^)1 

i =1,2,3 

VxVj 

= 2^7 [(* - x ))* + (y~ vj)$ + z ( z < + Az ~ z )] 

j = 4,5,6 ( 6 ) 

V x V* 

= 2 ^ e [(®*2 - **:l)z + (Vk 2 - Vkl)y] 

k = 7,8,9 


For FEM implementation, the following FEM elementary entries 


/>'„ = ^VxV.-VxV.ff (7) 

Q' m „ = J v V m V„dV ( 8 ) 


need to be explicitly expressed on the assumption of isotropic material prop- 
erty. To this end, we follow the notation defined in (5) and (6), where 
i, *'=1,2,3 represent the top triangle edges, j,j’= 1,2,3 denote the bottom tri- 
angle edges and fc, /c'=l,2,3 stand for the vertical three edges. It is found 
that (7) and (8) can be analytically evaluated and we tabulate the results as 
follows 
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Qn' 

Qkk‘ 

Qij 

Qik 


(Azf 


Cjj'Djji 


AzS e T k k> 

(Az ) 3 

Qa = ^j-CaDii 
Qki — Qjk = Qkj = 0 


(16) 

(H) 

(18) 

(19) 


where 
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1/6 fork = k'\ 1/12 fork^k' 
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The related quantities in the above list of the expressions are defined as 


s e = 

Is. dxdy 

sx = 

J x dxdy 

SY = 

/ y dxdy 
Js* 

SXX = 

J x 2 dxdy 

SYY = 

/ y 2 dxdy 
JS C 

SXY = 

/ xy dxdy 


and these integrals can be analytically evaluated and the results will be given 
in the appendix. 


Shape Functions: Type Two 

Edge-based shape functions for triangular elements can be derived and ex- 
pressed in various ways and we demonstrate here another type of shape func- 
tions. Either of these two sets may be employed for FEM implementation 
and the solutions from both FEM implementations are shown to be identical. 

We start from the nodal basis functions for 3-D triangular prism elements. 
It is well known that the nodal basis functions for the prism shown in figure 
2 may be simply expressed in terms of the simplex coordinates of triangles 
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Cl , C 2 , and ( 3 = I-C 1 -C 2 , where (. (i=l,2,3) are equal to the ratio where 
Si is the area of the triangle formed by an inner point and the two corners 
opposite to the i th corner. This definition provides a geometrical explanation 
as to why C is equal to 1 at node i and zero at the other two nodes. The 
difference of the functions for the top and the bottom triangle is taken into 
account by introducing a linear dependence in z so that the node shape 
functions associated with the upper triangle vanishes at the lower triangle. 


for nodes on top surface 

Af, = i-(z - *)C, 

N 2 = -T — (* - Zc) C2 

\ z 

* = Tz {z - 2t)C3 


for nodes on bottom surface 
N 4 = A (z c + A z- z)Ci 

Ns = 4~(Zc + A * - z )< 2 

T 

N 6 = -t ~{ z c + A z- z )( 3 
A z 


(21) 


The corresponding edge basis functions can be constructed using the 
above simplex coordinates. Specifically, the variations of functions along 
the z axis remain the same as in (21). However, in the x-y plane, we let 


S i = G,Vfo-&VCi„ i = l,2,3 (22) 

where *j, i 2 represent node 1 and node 2 of the i th edge, whose direction 
points from node 1 to node 2. There are several interesting properties as- 
sociated with this type of functions. First, since the simplex coordinates 
£ m , m = 1,2,3 are linear functions with respect to x and y, the gradient of 
the coordinates become constants. Specifically, 

VCm = z X (23) 


where l m is the length of the m th edge opposite to node m and S e is the 
triangle area, as usual. Also, it can be shown that the relation 


( 1 rG edge i 

1 0 r€ other edges 


(24) 


holds for r at any location of the edges, taking advantage of the features of 
the simplex coordinates. Clearly, (22) is also divergence-less, since V • S, = 0. 
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It is again straightforward to write down the edge-based shape functions for 
triangular prisms. They are given, by inspection, as 


V« — ^^(z 2c )(C<i VCi 2 Ci ’2 VCi, )j i— 1,2,3 

Vj = A~z (Zc + AZ ~ ~ 0.V0J, 3 = 4,5,6 

V fc = Kk k= 7,8,9 


(25) 


It is noted that \ k (k= 1,2,3) are the same as those associated with the type 
one shape functions. 

The curl of (25) is explicitly given by 


VxVi = AL {Cm z X VCi 2 - C. 2 z x VCm + 2 (z - z c )VC H x VC, 2 } (26) 

V X Vj = A~z {CnZ X V ° 2 " Cj2 ^ X VCjl + 2(Z ~ Zc)VCjl X V ^>( 27 ) 

V x Vjt = V( k x £ (28) 


and these expressions can be used to evaluate the matrix elements (7) and 
(8). The integral relations needed for this purpose are given in the appendix. 
Below we simply list the analytical expressions for the FEM entries: 
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Qik = Qh = Qjk — Qk] — 0 


( 39 ) 


12 ^ 

n n r 22 + n 22 r n — ni 2 T 2 i — n 2 j7i 2 
(|^ 2 (n„n 22 -n 12 n 21 ) 

VCi m -VCj n m,n = 1,2 

and 

/or = jn 

/or i m ± j n 

Also, the length of i th edge, for instance, is denoted by U whereas l k represents 
the length of the edge opposite to the vertex k. 

It can be seen that the two sets of shape functions have different geo- 
metrical meanings and both lead to closed form matrix entries. Also, the 
complexity of two formulations is about the same (careful observation shows 
that type-two is a bit more concise in form but type-one may be a bit easier 
to program) and the solutions rendered by their corresponding FEM imple- 
mentations are identical. 

III. Mesh Generation And Code Interface 

Triangular prismatic elements are practically attractive for conformal an- 
tenna simulations because they simplify mesh generation procedures rather 
substantially. Basically, once the surface triangular mesh of the printed an- 
tenna aperture is completed, the volume mesh can be generated automati- 
cally by growing the prismatic elements vertically into the cavity. Inhomo- 
geneities, layers and parasitic elements can be readily incorporated into the 
triangular prism mesh without difficulty. 



where 
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Typically, upon completion of the mesh, the following parameters/tables 
must be compiled for use by the finite element solver. They are: 

• node numbers and coordinates on top surface; 

• table of element numbers and corresponding nodes on top surface; 

• table for element and corresponding edge numbers and node numbers; 

• group of nodes at the outer boundary of top surface; 

• group of elements in the slot(s)/aperture(s); 

• additional group(s) of nodes (on surface) for interior PEC’s. 

Other information, such as number of layers along z axis, thickness and 
material properties (e r ,p r ) associated with each dielectric layer, as well as 
frequency scanning range, output information (scattering RCS/radiation pat- 
tern or input impedance) and iterative solution tolerance, should be inputted 
interactively. 


IV. Preliminary Results 

To validate our recently developed finite element-boundary integral code, we 
chose an annular ring slot antenna, whose tetrahedral mesh is available. Both 
RSC and radiation patterns are tested for the antenna. The configuration 
used for this testing is shown in figure 3, where the radius and the depth 
of the cavity are 12.35 cm and 3 cm, respectively; the ring has a radius 7.7 
cm and the slot width is 0.75 cm; the dielectric constant is 1.35. In both 
scattering and radiation cases, we chose the operational frequency to be 0.8 
GHz. 

Figure4 shows the comparisons of the RCS pattern with the computed 
results using the tetrahedral FE-BI code, whereas figure 5 demonstrates both 
the co-polarization and cross-polarization radiation patterns. They all are 
in good agreement. The mesh generation procedure for the prismatic element 
implementation is significantly simplified. 

Our near future work is to test various spiral antennas using this code and 
and model the input impedance for different feed structures. Also, the impact 
of the dielectric substrate/superstrate on the performance and geometry of 
antennas, such as polarization control and structural miniaturization will be 
especially focused. 
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a = 12.35 cm 



b = 0.75 cm 
p 0 = 7.7 cm 
0.7 < f < 1 GHz 


d = 3 c m| [ "[ 


Figure 3: The geometry of the annular ring slot antenna used for comparisons 
of the prismatic FE-BI code with the tetra FE-BI code. 

Appendix 

Certain integral relations have been employed in the report for analytical 
evaluation of the elementary FEM entries. These are often found in FEM 
books but they are usually in local coordinates. Since the expressions in 
global coordinates are also needed, we tabulate them in this appendix without 
proof. 

Assume the three nodes i,j and m of a triangle are in counterclockwise 
rotation, their global coordinates are given by (A,, VJ-), (Xj, Yj), (X m , Y m ). we 
have, 


S e = [ dxdy = ^ 

1 

Xi 

y i \ 

1 

Xj 

y j I 

ys<= 2 

1 


Vm ) 


SX = f x dxdy — — (Xi + Xj + X m ) 

JS e o 

r S e 

SY = / ydxdy = —(Y, + Y : + Y m ) 

JS c o 
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Figure 4: Comparisons of the bistatic (co-pol) RCS patterns computed using 
the tetra FE-BI code and the prismatic FE-BI code. The normal incident 
plane wave is polarized along the 0=0 plane and the observation cut is 
perpendicular to that plane. 



Figure 5: Illustration of x-pol and co-pol radiation patterns from the ring 
slot antenna shown in figure 3. The solid lines are computed using the 
tetrahedral FE-BI code whereas the dotted lines are computed using the 
prismatic FE-BI code. The excitation probe is place at the point (y=0) in 
the slot marked in figure 3. 
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SXX = X 2 dxdy = y^- {(A; + A'j + A m ) 2 + (A ' 2 + A 2 + A^) } 

SET = J s ' y 2 dxdy = ^ {( K + Y- + Y m ) 2 + (V ' 2 + Y/ + >„) } 

SXY = f xydxdy = ^r{(A, + X : + A„,) (K + Yj + Y m ) 

Js* 

+ (XiYi + A jYj + A m Y m )} 


Also, if the simplex coordinates for these three nodes are denoted by 
(C»>CnCm)i then 
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j sc cr c] Cm dxd v 
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where p, q, r are all non-negative integer numbers. Since the elements we used 
for the FEM implementation are linear, the associated integrands involve only 
linear or quadratic functions of C- Therefore, 
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Abstract 

A new artificial material absorber for truncating finite element meshes is 
investigated. The interface of the absorber is made reflectionless by choosing e r 
and fi T to be complex diagonal tensors. With some loss, a metal backed thin ab- 
sorber layer is then sufficient for terminating the mesh. This scheme is simpler 
to implement than conventional absorbing boundary conditions and offers the 
potential for higher accuracy. In this paper, we investigate the effectiveness of 
this anisotropic absorber on the basis of results obtained for problems in prop- 
agation (waveguide and printed circuits) and scattering. Based on the results, 
specific recommendations are given for selecting the absorber loss parameters 
to optimize its performance for a given thickness. 


I. Introduction 

One of the most important aspects of finite difference and finite element implemen- 
tations is the truncation of the computational volume. An ideal truncation scheme 
must ensure that outgoing waves are not reflected backwards at the mesh termina- 
tion surface, i.e. the mesh truncation scheme must simulate a surface which actually 
does not exist. To date, a variety of non-reflecting or absorbing boundary conditions 
(ABCs) have been employed for truncating the computational volume at some dis- 
tance from the radiating or scattering surface, and applications to microwave circuits 
and devices have also been reported. The ABCs are typically second or higher order 
boundary conditions and are applied at the mesh termination surface to truncate the 
computational volume as required by any PDE solution. Among them, a class of 
ABCs is based on the one-way wave equation method [1, 2] and another is derived 
starting with the Wilcox Expansion [3, 4]. Also, higher order ABCs using Higdon’s 
[5, 6]formulation and problem specific numerical ABCs have been successfully used, 

1 EM CAD Lab, ECE Department, Worcester Polytechnic Institute, Worcester, MA 01609 
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particularly for truncating meshes in guided structures [7]. There are several diffi- 
culties with traditional ABCs. Among them are accuracy limitations and a need to 
enforce them at some distance from the perturbance (leading to a large computational 
volume). Also, although higher order ABCs are considered more accurate, they com- 
promise the convergence of the solution and, furthermore, the resulting system is not 
as easily amenable to parallelization. 

An alternative to traditional ABCs is to employ an artificial absorber for mesh 
truncation. Basically, instead of an ABC, a thin layer of absorbing material is used to 
truncate the mesh, and the performance for a variety of such absorbers has been con- 
sidered [8], [9]. Nevertheless, these lossy artificial absorbers (homogeneous or not) still 
exhibit a non-zero reflection at incidence angles away from normal. Recently, though, 
Berenger [10] introduced a new approach for modeling an artificial absorber that is 
reflectionless at its interface for all incidence angles. In two dimensions, his approach 
requires the splitting of the field components involving normal (to the boundary) 
derivatives and assigning to each component a different conductivity. In this manner 
an additional degree of freedom is introduced that is chosen to simulate a reflection- 
less medium at all incidence angles. Provided the medium is lossy, this property is 
maintained for a finite thickness layer. Berenger refers to the latter as a perfectly 
matched layer(PML) and generalizations of his idea to three dimensions have already 
been considered [11], [12]. Also, implementations of the absorber for truncating finite 
difference-time domain(FDTD) solutions have so far been found highly successful. 
Nevertheless, it should be noted that Berenger’s PML does not satisfy Maxwell’s 
equations and cannot be easily implemented in finite element (FEM) solution. 

In this paper we examine a new anisotropic(uniaxial) artificial absorber [13] for 
truncating FEM meshes. This artificial absorber is also reflectionless at all incidence 
angles. Basically, by making appropriate choices for the constitutive parameter ten- 
sors, the medium impedance can be made independent of frequency, polarization, 
and wave incidence angle. A PML layer can then be constructed by introducing suf 
ficient loss in the material properties. The implementation of this artificial absorber 
for truncating finite element meshes is straightforward and, moreover, the absorber 
is Maxwellian. Below, we begin with a brief presentation of the proposed artificial 
absorber, and this is followed by an examination of the absorber’s performance in ter- 
minating guided structures and volume meshes in scattering problems. Results are 
presented which show the absorber’s performance as a function of thickness/frequency 
and for different loss factors. 


II. Formulation 

Consider the waveguide, shielded microstrip and scatterer shown in Figures 1 and 
2. We are interested in modeling the wave propagation in these structures using 
the finite element method. For a general anisotropic medium, the functional to be 
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minimized is 


?= j V x E ■ (£-' • V x E) - k% • E - El iV 

-jL. E x (£' • V x E) • dS, (1) 

in which Ji T and c r denote the permeability and permittivity tensors whereas E is 
the total electric field in the medium. The surface integrals over S, n and S out must 
be evaluated by introducing an independent boundary condition and the ABC serves 
for this purpose but alternatively an absorbing layer may be used. An approach to 
evaluate the performance of an absorbing layer for terminating the FE mesh is to 
extract the reflection coefficient computed in the presence of the absorbing layer used 
to terminate the computational domain. In this study we consider the performance of 
a thin uniaxial layer for terminating the FE mesh for certain wave guided structures 
and scatterers. Such a uniaxial layer was proposed by Sacks et.al. [13] who considered 
the plane wave reflection from an anisotropic interface (see Figure 3 ). If and f r 
are the relative constitutive parameter tensors of the form 

/ a 2 ° ° \ 

Mr = e r = 0 b 2 0 (2) 

\ 0 0 c 2 J 

the TE and TM reflection coefficients at the interface (assuming free space as the 
background material) are 

cosOi — <J~^cos6 t 
cosOi + \f^-cos6 t 

V “2 ZO\ 

yJ^COs6 t — COs6i V 

cosOi + \p^cosB t 

and by choosing a 2 = b 2 and c 2 = ^ it follows that R TE = R™ = 0 for all incidence 
angles, implying a perfectly matched material interface. If we set a 2 = a — jfl, the 
reflected field for such a metal-backed uniaxial layer is 

\R(ei)\ = e ~ 2(3ktcos e ' ( 4 ) 

where t is the thickness of the layer and 0, is the plane wave incidence angle. The 
parameter ak is simply the wavenumber in the absorber. Basically, the proposed 
metal backed uniaxial layer has a reflectivity of -30 dB if /3tcos$i = 0.275A or -55dB if 
/3tcos9i = 0.5A, where A is the wavelength of the background material. The reflection 
coefficient (4) can be reduced further by backing the layer with an ABC rather than 
a PEC. However, the PEC backing is more attractive because it eliminates altogether 
the integrals over the surfaces. Clearly, although the interface is reflectionless, the 
finite thickness layer is not, and this is also true for Berenger’s PML absorber. 


R TE = 

rtm = 
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Below we present a number of results which show the performance of the proposed 
uniaxial absorbing layer a s a function of the parameter /?, the layer thickness t and 
frequency for the guided and scattering structures shown in Figures 1 and 2. We 
remark that for the microstrip line example it is necessary to set a 2 = e r f>(c* - j(3) for 
the permittivity tensor and a 2 = Hrb{ Q ~ jft) f° r the permeability tensor, where e r t 
and fi r b are the relative constitutive parameters of the background material (i.e. the 
substrate). 


Results 

Rectangular Waveguide 

Let us first consider the rectangular waveguide shown in Figure 1. The guide s cross- 
section has dimensions 4.755 cm x 2.215 cm and is chosen to propagate only the T E 10 
mode. It is excited by an electric probe at the left, and Figure 4 shows the mode 
field strength inside the waveguide which has been terminated by a perfectly matched 
uniaxial layer. As expected, the field decay inside the absorber is exponential and 
for (3 values less than unity the wave does not have sufficient decay to suppress 
reflections from the metal backing of this 5cm layer. Consequently, a VSWR of about 
1.1 is observed for /3 = 0.5. However, as /? is increased to unity, the VSWR is nearly 
1.0 and the wave decay is precisely given by e - ^ 009 **, where d is the wave travel 
distance measured from the absorber interface and here 0,- = 44.5°. Nevertheless, 
when (3 is increased to larger values, the rapid decay is seen to cause unacceptable 
VSWR’s. One is therefore prompted to look for an optimum decay factor f3 for a 
given absorber thickness and Figure 5 provides a plot of the TE\q mode reflection 
coefficient as a function of (3. Figure 5 is typical of the absorber performance and 
demonstrates its broadband nature and the existence of an optimum value of (3 for 
minimizing the reflection coefficient. Basically, the results suggest that (3 must be 
chosen to provide the slowest decay without causing reflections from the absorber 
backing. That is, the lowest reflection coefficient is achieved when the entire absorber 
width is used to reduce the wave amplitude before it reaches the absorber’s backing. 
As expected, this optimum value of (3 changes with frequency but the broadband 
properties of the absorber are still maintained since acceptable low reflections can 
still be achieved for unoptimized (3 values. For example, in the case of / = 4.5 Ghz 
the optimum value of /3 = 1 gives a reflection coefficient of -45dB whereas the value 
of (3 = 3 gives a reflection coefficient of -37 dB which is still acceptable. It should be 
noted though that setting /? = 3 allows use of an absorber which is about 2cm or 0.3 
free space wavelengths. 

Not surprisingly (see (4)), for this example, the value of a does not play an im- 
portant role in the performance of the absorber and this is demonstrated in Figure 6. 
As seen, setting a = (3 gives the same performance as the case of a = 1 shown in 
Figure 5. Our tests also show that other choices of or give the same absorber perfor- 
mance. However, it is expected that a will play a role in the presence of attenuating 
modes and it is therefore recommended to choose a — (3 to ensure that all modes are 
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absorbed. 

Microstrip line 

The performance of the perfectly matched uniaxial layer in absorbing the shielded 
microstrip line mode is illustrated in Figure 7 where the reflection coefficient is plotted 
as a function of 0. In this case, the microstrip line is extended into the 5cm absorbing 
layer up to the metallic wall. Similarly to the waveguide, we again observe that an 
optimum / 3 value exists and it was verified that (in the absorber) the quasi-TEM mode 
has the same attenuation behavior as shown in Figure 4. The reflection coefficient at 
the optimum 0 = 1 is now —28.5 dB and if better performance is required, a thicker 
absorbing layer may be required. Again, as in the case of the waveguide example, 
the value of a is of small importance to the performance of the absorber and this is 
demonstrated in Figure 8. 

Band Eliminator 

Another guided structure that was considered is the stripline band eliminator shown in 
Figure 9. The band eliminator consists of a circular disk 1.64 cm in radius, placed on 
a substrate having a relative dielectric constant of 2.4. The disk is directly coupled 
to 50 Cl striplines and to take advantage of symmetry, only half of the problem’s 
domain is used in the FEM computation. The thickness of the absorber used in this 
calculation is 1.5cm, and in Figure 10 we present the insertion loss results from 2GHz 
to 4GHz. The values of a and 0 for the absorber were both set to 3 and, as seen, our 
calculations compare very well with measured data and reference calculations given 
in [14], thus, demonstrating the accuracy of the absorber termination. 

Sphere Scattering 

As another example we examined the scattering by a metallic sphere illuminated by a 
plane wave (see Figure 2). The unknown quantity in the FEM implementation was the 
electric field, and the excitation due to the incident plane wave was incorporated by 
enforcing the boundary condition E® = — E[ on the surface of the sphere. The FEM 
mesh was terminated by a cubical 0.1 A 0 thick metal backed absorber whose interface 
was placed only 0.05A o away from the sphere’s surface. Based on the attenuation 
curves shown in Figure 4, a value of 0 = 4 or slightly greater is needed to suppress 
reflections from the absorber’s backing. Indeed, lower values of 0 did not perform well 
for this application. The bistatic scattering pattern computed with a mesh density of 
14 edges per wavelength is shown in Figure 11 and is within ldB of the exact result. 
Better agreement would clearly require thicker absorber and a lower 0 value. Again, 
the value of a does not have appreciable effect on the absorber’s performance. About 
56,000 unknowns were used for this simulation and the system converged in about 
2000 iterations. We also remark that higher sampling densities did not appreciably 
improve the accuracy of the result or the relative convergence rate. 


Conclusions 

A new method of mesh truncation based on an artificial absorbing layer with anisotropic 
material properties was implemented for terminating the FEM mesh. Explicit recom- 
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mendations were given for choosing the absorbers thickness and attenuation properties 
to optimize performance. It was demonstrated that good absorption can be achieved 
with relatively thin absorbing layers (a small fraction of a wavelength) placed close 
to the computational domain. Also, because there is no surface integral or ABC at 
the mesh truncation boundary, the method is easily implemented in FEM solutions 
and is well suited for parallel computations. 
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(a), waveguide 
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h ■ 021 cm 
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(b). Stripline 

Figure 1: A rectangular waveguide (a) and a microstrip line (b) truncated using the 
perfectly matched uniaxial absorbing layer. 



Figure 2: Geometry of sphere scattering example. 
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Figure 3: Plane wave incidence on an interface between two diagonally anisotropic 
half- spaces. 



Figure 4: Field values of the TEi 0 mode inside a waveguide terminated by a per- 
fectly matched uniaxial layer. The absorber is 10 elements thick from the 71 4t to the 
80 t/l element and each element was 0.5 cm which translates to about 13 samples per 
wavelength at 4.5 GHz. 


50 





Figure 5: Reflection coefficient vs /? (a = 1) for the perfectly matched uniaxial layer 
used to terminate the waveguide shown in Figure 1(a). 



Figure 6: Reflection coefficient vs /?, with a = f3, for the perfectly matched uniaxial 
layer used to terminate the waveguide shown in Figure 1(a). 
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Figure 7: Reflection coefficient vs /? with a= 1, for the shielded microstrip line termi- 
nated by the perfectly matched uniaxial layer. 



Figure 8: Reflection coefficient vs (3 with a = (3, for the shielded microstrip line 
terminated by the perfectly matched uniaxial layer. 
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Figure 10: Insertion loss for the band eliminator shown in Figure 9. 
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Frequency (GHz) 



Bistatic RCS (14 elements per wave length) 



Figure 11: Bistatic scattering for a 0.4A diameter sphere for different values of (3. The 
FEM mesh density was approximately 14 elements per A and was terminated by a 
cubical thin absorber whose inner and outer diameter was 0.5A and 0.7A, respectively. 
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Abstract 

We examine the performance and optimal selection of a new uniaxial artificial 
absorber layer for truncating finite element meshes. The absorber is reflectionless at 
its interface for all incidence angles but in practice a finite layer must be used, and 
as a result the absorber is no longer reflectionless. In this Letter we examine the 
absorber’s optimal selection of its thickness and dielectric properties in the context 
of a finite element implementation for a shielded microstrip line. 


Introduction 

A most important aspect of a finite element implementation is the termination of the 
computational domain. This is carried out by either using an absorbing boundary con- 
dition (ABC) [1, 2, 4, 3] or a thin layer absorber [5, 6]. However, the ABCs and thin 
layer absorbers considered so far have non— zero reflections for waves impinging at oblique 
incidences and this limits their accuracy and leads to requirements for larger computa- 
tional domains. Recently, Sack et al. [7] proposed a new anisotropic (uniaxial) absorber 
for truncating finite element meshes. This absorber is indeed reflectionless (i.e. perfectly 
matched at its interface) for all incidence angles. However, in practice a finite thickness 
absorber must be employed, as shown in figure 1. As a result, although the absorber has 
a perfectly matched interface (PMI), the layer is no longer reflectionless due to reflection 
from the back of the layer. It is therefore of interest to optimize the layer’s absorptivity 
in terms of its thickness and material properties. Of interest is to also examine the ef- 
ficiency of the absorber in actual numerical simulations, since to date no application of 
the absorber’s performance and optimization has yet been published. In the following, we 
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examine the absorber’s performance and optimization when used to terminate the finite 
element mesh of a shielded microstrip line. 


Uniaxial Absorbing Layer 


Consider the uniaxial absorbing layer shown in Figure 1 whose relative permittivity and 
permeability are given by 



( a 

0 

°\ 

( a 
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°\ 

Hr = Vrb 
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( 1 ) 


where (e r (,,^ r f>) refer to the relative permittivity and permeability of the background 
medium. The corresponding TE and TM reflection coefficients at the interface are given 
by 


R TE 

R™ 


cosd, — yJ^COsOt 
COS0i + yJ^COsOt 
yJ^COsOt — COS&i 
COS0i + yj\cos0 t 


( 2 ) 


and by choosing a = b and c = 1/a, it follows that R TE = R™ = 0 for all incidence angles. 
This implies a perfectly matched interface (PMI) and the transmitted wave decays as 
e -pkdcos6, , w } iere j i s the distance from the interface, k = u>y/fi 0 toy/fJ, r ^r is the wavenumber 
of the background material and we have set a = a — jft. Clearly, unless the transmitted 
wave is completely absorbed while traveling in the metal-backed absorber, the layer will 
have a nonzero reflected field. It is therefore of interest to optimize the absorbing layer 
performance through an appropriate selection of its parameters. 


Application to Shielded Microstrip Line 

Consider the microstrip line shown in Fig. 2. The line is excited by a vertical probe at its 
left and is terminated to the right by a thin PMI layer of the uniaxial absorber. The length 
of the line is d+t=12cm and a finite element implementation with the brick element [8] 
was used for modeling the shielded structure. We also note that for this implementation 
the microstrip line was extended into the 5cm absorbing layer. Figure 2 shows the vertical 
field magnitude below the microstrip line for different ft values and a = 1. Only the fields 
across and near the absorbing layer are depicted since our interest is to examine the field 
behavior in that region. As expected the field decay inside the absorber is exponential and 
for ft less than unity the wave does not have sufficient decay to suppress reflections from 
the metal backing of the 5cm PMI layer. As ft increases to unity, the VSWR decreases 
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to about 1.05 but for (3 values above 3 we again observe deterioration of the VSWR. One 
is therefore prompted to seek an optimum decay factor /3 for a given absorber thickness. 
Our study led to the conclusion that (3 should be selected so that the entire absorber 
thickness is utilized in attenuating the field to less than 5 percent of its original value. 
For the chosen 5cm PMI layer, a value of (3 = 1 or so is appropriate at a frequency of 
4.5GHz for a 30 dB reduction. Lower reflectivities can be achieved by increasing the layer 
thickness and changing (3 appropriately. It turns out that a plays very little role in the 
performance of the absorber since the curves in Figure 2 remain nearly unchanged as a 
is increased above unity. For large (3 values the decay becomes very rapid and the PMI 
interface begins to resemble a magnetic wall. This behavior is due to the numerical FEM 
implementation and stems from the inaccurate modeling of the rapidly decaying fields in 
the layer. It is basically more efficient to use a thicker absorber layer with a lower f3 value 
and use nominal sampling (10 to 15 per wavelength) rather than a thinner absorber with 
high (3 values and consequently excessive sampling requirements. Our study shows that 
the absorber fields with 1 < f3 < 2 can be accurately modeled using cell sizes on the order 
of 0.05 guided wavelengths. 

The accuracy of the PMI layer for circuit parameter computations is illustrated in 
Figure 3. It is seen that use of the optimized 5cm PMI layer, with a = 1 and (3=1, 
yields very accurate input impedance values. The shown microstrip line impedances were 
computed by measuring the vertical field at the probe’s location without a need to extract 
the VSWR which is often difficult with unstructured finite element meshes. A recent 
measurement technique has also permitted the validation of the computed tangential 
fields (to the microstrip) above the dielectric. As seen from Figure 4, the calculated and 
measured fields are nearly overlaid and this is of importance since inexact mesh truncation 
schemes are prone to inaccuracies in the near zone. Note that the shielded microstrip line 
dimensions for the data in Fig. 4 are different from those corresponding to Fig. 3. 


Conclusion 

We examined the performance and optimal use of a new uniaxial absorber layer with a 
perfectly matched interface(PMI) for terminating finite element meshes. It was found that 
the loss parameter of the PMI layer must be selected to have a nominal value and make 
full use of the entire layer thickness for field decay. Moreover, it was demonstrated that 
such a PMI layer is extremely accurate for truncating finite element meshes and rivals the 
Berenger’s PML [10] layer which was recently introduced for truncating finite difference 
-time domain meshes. It is therefore expected that the PMI layer will find wide-spread 
use for truncating finite element meshes. 


58 



References 


[1] L.Halpern and L.N. Trefethen, “Wide-angle one-way wave equations,” J. Acoust. 
Soc. Amer. Vol. 84, pp. 1397-1404, 1988 

[2] R.L. Higdon, “Absorbing boundary conditions for acoustic and elastic waves in strat- 
ified media,” J. Comp. Phys. Vol. 101, pp. 386-418, 1992 

[3] J.Fang, “ABCs applied to model wave propagation in Microwave integrated- cir- 
cuits,” IEEE Trans. MTT, Vol. 42, No. 8, pp. 1506-1513, Aug. 1994 

[4] A. Chatterjee and J.L.Volakis, “Conformal absorbing boundary conditions for the 
vector wave equation,” Microwave Opt. Tech. Lett. Vol. 6, pp. 886-889,1993 

[5] T. Ozdemir and J.L.Volakis, “A comparative study of an absorbing boundary condi- 
tion and an artificial absorber for truncating finite element meshes,” Radio Science 
Vol. 29, No. 5, pp. 1255-1263, Sept. 1994 

[6] C. Rappaport and L. Bahrmasel, ”An absorbing boundary condition based on ane- 
choic absorber for EM scattering computation,” J. Electromagn. Waves Appl., Vol. 
6, No. 12, pp. 1621-1634, Dec. 1992. 

[7] Z.S. Sacks, D.M. Kingsland, R. Lee and J.F. Lee, “A perfectly matched anisotropic 
absorber for use as an absorbing boundary condition,” submitted for publication 

[8] J.L. Volakis, J.Gong and A. Alexanian, “Electromagnetic Scattering From Microstrip 
Patch Antennas and Spirals Residing in a Cavity,” Electromagnetics, Vol. 14, No. 1, 
pp 63-85, 1994 

[9] B.C. Wadell, “Tramsmission Line Design Handbook,” Artech House, Boston, 1991 

[10] J.P. Berenger “A Perfectly Matched Layer for the Absorption of Electromagnetic 
Waves,” J. Comp. Physics , Vol. 114, pp. 185-200, 1994. 


59 




60 





Figure 2: 



10 100 1000 
Position (microns) 

Figure 3: Comparison of the near field as a function of locations above the microstrip 
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1 Abstract 

This paper considers the hybridization of the finite element and high 
frequency methods for predicting the radiation pattern of printed 
antennas mounted on aircraft platforms. The finite element method 
is used to model the cavity-backed antennas whereas the interac- 
tions between the radiators and the substructures are treated via a 
high frequency technique such as the GTD, PO/PTD or SBR. We 
present comparisons between measurements and calculations along 
with a qualitative description of the employed finite element and 
high frequency codes. 
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2 Introduction 


Over the past few years the finite element method (FEM) has been 
applied to various scattering [l]-[3] and antenna [4]-[6] analysis prob- 
lems with remarkable success. One of the main advantages of the Fi- 
nite Element Method (FEM) is its inherent geometrical adaptability 
(see Figure 1) and ease in handling materials. Consequently rather 
complex configurations can be readily handled without a need for 
code rewriting or modification. Typically, the FEM solver is inter- 
faced with sophisticated mesh generators [7]-[8], which can provide 
the necessary volume mesh using a variety of tessellation elements 
such as bricks, tetrahedrals or prisms (see Figure 2). This permits 
an accurate geometrical modeling of the structure which is essential 
for antenna related problems since fine geometrical features play an 
important role in the characterization of the antenna parameters. 
However, as is often the case with robust numerical simulations, the 
combined modeling of large and small geometrical features within a 
single mesh or computational domain leads to inefficiencies and un- 
acceptable CPU requirements. An example of this situation is the 
problem of antenna radiation in the presence of a large platform such 
as an aircraft, a satellite or other airborne and land vehicles. In this 
case, the antenna requires accurate modeling of its fine geometrical 
features and material parameters (see Figure 3) using, for example, 
the FEM. However, modeling of the large aircraft platform using 
the FEM requires several million degrees of freedom, thus, making 
it unattractive for such an implementation. Instead, it is more ap- 
propriate to model the vehicle substructure using a high frequency 
method such as the Geometrical Theory of Diffraction (GTD) [9] or 
the Physical Optics(PO)/Physical Theory of Diffraction (PTD) [10] 
and Shooting and Bouncing Ray(SBR) method [11]. 

To overcome the inefficiencies of a single formulation, in this paper 
we present a simple hybridization of the finite element and high fre- 
quency methods. Basically, the finite element method is employed 
for the analysis of the antenna in the absence of substructures such 
as wings, fins, engines, or surface details. The pattern or antenna 
aperture currents generated by the FEM code are then used as the 
source(s) of the high frequency code. Two high frequency codes 
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were investigated in examining the effectiveness of this hybridiza- 
tion. One of the codes is based on the GTD and the other high 
frequency code combines the PO/PTD and SBR to compute the 
interactions of the antenna pattern with the surrounding substruc- 
tures. Basically, the proposed hybridization takes advantage of the 
robust FEM formulation to obtain an accurate characterization of 
the antenna in isolation whereas the high frequency technique is used 
for modeling the interactions with the large features of the substruc- 
ture^ Figure 4). The effectiveness of this hybridization is eval- 
uated by comparing the results of the analysis with measurements 
corresponding to a patch antenna radiating on a finite cylinder in 
the presence of a flat plate attached to the cylinder as illustrated 
in Figure 5. In the following sections we first give a qualitative de- 
scription of the finite element and high frequency codes, and then 
proceed with a more detailed description of the hybridization proce- 
dure followed by a comparison of measurements with calculations. 

3 Finite Element Code Description 

Several finite element codes have been developed at the University 
of Michigan for the analysis and design of printed antenna configu- 
rations. Typically, the printed antenna configuration is assumed to 
be recessed in some metallic or coated platform and the codes dif- 
fer in the element used for the tessellation of the antenna, the type 
of platform assumed in the analysis (planar, cylindrical or doubly 
curved) and the closure condition employed for terminating the fi- 
nite element mesh as illustrated in Figure 6. The following codes 
are available for antenna radiation and scattering analysis. 

FEMA-CYL: This code is specialized to cavity-backed antennas 
recessed in a metallic cylindrical platform(see Figure 6). The FEM 
is used to model the cavity region and the boundary integral or ab- 
sorbing boundary conditions (ABC) are employed for truncating the 
mesh. In order to simplify the user interface, cylindrical shell ele- 
ments are used for the discretization of the cavity and the aperture 
(see Figure 2). The result is a structured mesh and the user inter- 
face is reduced to specifying the dimensions of the antenna. Multiple 
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patches and patch arrays (consisting of individual cavity-backed ele- 
ments or elements on a single substrate) can be considered but their 
geometry is restricted by the generated grid. This situation is similar 
to finite difference codes, where the model’s boundary is modified to 
fit the geometry. Consequently, although this code is user-friendly, 
it is not appropriate for computing the radiation parameters of cir- 
cular patches or spirals. However, scattering parameters are not as 
sensitive to minor geometrical modifications and therefore FEMA- 
CYL can still be used for computing scattering from non-rectangular 
printed antenna configurations and cavities. 

The theoretical background of FEMA-CYL is described in [12] and 
[13], and the code has been validated for antenna and scattering ap- 
plications. An example calculation using FEMA-CYL is illustrated 
in Figure 7. The measured data given in Figure 7 were obtained 
at NASA-Langley [14] with the patch antenna placed on the cylin- 
drical section of the ogive+cylinder structure as illustrated in the 
Figure. For this calculation, the cylinder Green’s function of the 
second kind was employed for truncating the finite element mesh to 
generate a combined finite element- boundary integral system. The 
latter is solved using the biconjugate gradient [15] or QMR [16], [17] 
iterative solver and the FFT is used to speed-up the matrix- vector 
product evaluation in the solver [18]. Consequently, in spite of the 
fact that the boundary integral matrix is fully populated, the code’s 
computational and memory requirements remain at O(N). The sys- 
tem solution generates the fields within and on the surface of the 
antenna cavity. For the radiation and scattering pattern calcula- 
tions only the surface electric fields are used to generate equivalent 
magnetic currents. These are subsequently substituted in the ra- 
diation integral with the appropriate platform Green’s function to 
generate the antenna parameters. For antenna excitation, a probe 
or a coaxial cable feed model is employed and in the case of scatter- 
ing the excitation becomes the aperture fields due to the incoming 
field. When a probe feed model is used, the input impedance is 
the line integral of the electric field along the length of the probe 
divided by the magnitude of the probe current. When a coaxial 
feed is employed, the excitation is introduced by setting the fields 
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coinciding with the edge elements bordering the opening of the coax 
cable equal to AV / A L, where AV is the given potential between 
the outer and inner surface of the cable and A L = b — a, where 
b and a correspond to the outer and inner radius of the cable, re- 
spectively. The center conductor is modeled by setting to zero the 
field unknowns that are associated with the edges coinciding with 
the conductor. The improved results of this feed modeling approach 
are shown in Figure 8 and the details of the modeling are described 
in [19]. 

FEMA-CYL can also be used to include the effect of superstrate ma- 
terials (see Figure 6) which may extend over the antenna platform. 
To avoid use of complicated Green’s functions, in this case the mesh 
is extended a fraction of a wavelength over the cavity’s aperture 
and an absorbing boundary condition (ABC) is employed for trun- 
cating the FEM mesh [10]. With this type of mesh truncation, the 
entire system is sparse and thus the memory and CPU requirement 
are again O(N) without a need to make use of the FFT. We have 
found that placement of the mesh at a distance of 0.3 wavelengths 
away from the coatings surface is sufficient to obtain the radiation 
pattern with reasonable accuracy as illustrated in Figure 7(b). The 
pattern is computed by introducing equivalent electric and magnetic 
currents placed at the surface of the dielectric coating and then prop- 
agating them using the free-space Green’s function. 

FEMA-TETRA: This code employs the same FEM formulation 
as FEMA-CYL for computing the parameters of printed antennas 
recessed in a metallic platform. In contrast to FEMA-CYL, this 
code employs tetrahedral elements for modeling the radiating struc- 
ture and bricks or shells for modeling the feed (see Figure 9). As 
a result, it incorporates maximum geometrical adaptability and has 
already been employed for the analysis of antenna+feed configura- 
tions beyond the capabilities of moment method codes. The FFT 
can also be invoked for carrying out the matrix- vector products as- 
sociated with the boundary integral portion of the system in the 
iterative solver, thus maintaining an O(N) memory requirement for 
the entire system as described in [4]. Its accuracy has already been 
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demonstrated for a vaxiety of antennas, including stacked patches 
and finite arrays in ground planes, printed and slot spirals and aper- 
ture antennas. In the case of aperture fed patch configurations such 
as that shown in Figure 10, the mesh generation is divided into the 
regions above and below the feeding slot. The mesh in each re- 
gion can employ different elements and this is important in avoiding 
meshing bottlenecks due to the narrowness of the slot. Two mesh- 
ing regions are then mathematically connected by enforcing electric 
field continuity across the slot. However, since the mesh elements in 
each regions are different, this condition is implemented by first re- 
lating the fields in the slot to the potential existing across the slot. 
Using this procedure, there is no need to have coincidence of the 
edges bordering the slot and thus different elements can be used in 
the regions above and below the slot aperture. The FEMA-TETRA 
code is interfaced with commercial solid- modeling packages such as 
SDRC I-DEAS and PATRAN [8] which are used to enter the geom- 
etry and generate the volume mesh. Using such geometry modeling 
packages and owed to the adaptability of tetrahedrals, nearly any 
printed antenna configuration can be modeled and analyzed using 
FEMA-TETRA. As can be expected, because FEMA-TETRA is not 
specialized to any antenna configuration all geometrical data must 
be entered through the solid- modeling package and this can be a 
time consuming task for inexperienced users. 

FEMA-PRISM: To avoid the time consuming task of volume 
mesh generation without compromising geometrical adaptability, 
this code provides an attractive compromise by making use of prisms 
(see Figure 2). Prisms have triangular faces at their top and bot- 
tom faces and can thus be used to model any patch configura- 
tion with sufficient geometrical fidelity. Typically, since the sub- 
strate/superstrate is of constant thickness, once the surface grid is 
constructed, the prismatic volume mesh can be generated by grow- 
ing the elements above and below the surface of the printed antenna 
as illustrated in Figure 11. Thus using prisms reduces the mesh 
generation process from a three-dimensional to a two-dimensional 
task. That is, the user needs to only generate the surface mesh in 
the aperture containing the printed antenna and the volume mesh 
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is then built automatically by specifying the substrate/superstrate 
thickness and information relating to the closure condition. FEMA- 
PRISM makes use of an ABC or metal-backed absorber for trun- 
cating the mesh at some distance from the aperture surface. The 
accuracy of this approach is illustrated in Figure 12 where the input 
impedance of a circular patch computed using a metal-backed ab- 
sorber closure is shown to be in good agreement with results based 
on the boundary integral closure. As in the other codes, the ra- 
diation and scattering patterns are computed by propagating the 
equivalent magnetic currents (existing at the cavity aperture) using 
the free-space Green’s function. Also, since a non-integral closure 
condition is used for truncating the mesh, no restriction is placed 
upon the curvature of the platform. Results generated by using the 
FEM codes refered to above will not, of course, include effects due to 
interactions between the antenna and the substructures residing on 
the same platform. To include these, the FEM codes are interfaced 
with the high frequency analysis packages described next. 

4 High Frequency Code Descriptions 

GTD/UTD Code: The hybridization of GTD and low frequency 
codes such as the moment method were considered in the early 1980s 
[20], [21] for scattering and antenna radiation analysis. Typically, the 
GTD was employed to account for diffraction contributions from 
edges and corners near the radiating elements and a review of some 
of this work was given more recently in [22]. In the case of scattering, 
hybrid GTD and moment method codes have been used to model 
small details on larger structures. Specific applications include uses 
of the moment method to model a small wire or a crack in the pres- 
ence of a large platform which is associated with several scattering 
centers caused by geometrical optics reflections and diffractions from 
edges and corners. For the most part these hybrid codes have been 
restricted to specialized implementations where the scattering sub- 
structure is hard-wired into the code and is characterized by a few 
scattering centers. Also, for radiation analysis available GTD codes 
have been restricted to include a small class of antenna elements 
such as narrow slots, monopoles and infinitesimal current elements. 


68 



A well developed UTD code, referred to as NEC-BSC was developed 
by Marhefka and Burnside [23] at the Ohio State University and was 
used in conjunction with the FEM codes for printed antenna pattern 
calculations on simple aircraft platforms. The UTD code computes 
the interactions between current elements and the aircraft substruc- 
ture using ray diffraction methods. In the employed UTD code the 
fuselage is modeled as a truncated cylinder or an ellipsoid and the 
attached wing and fins are modeled as combination of flat plates. 
The overall radiation pattern of the specified current element is then 
computed by adding the direct, geometrical optics, and diffraction 
contribution from all possible primary and secondary sources which 
arrive at the observer. Shadowing is also performed as part of the 
pattern calculation. In our implementation, the UTD code is sup- 
plied with the location and amplitude of the magnetic surface cur- 
rent elements as computed by the FEM code in the absence of the 
fuselage appendages. Consequently, these currents include all cou- 
pling effects internal to the antenna cavity and among the antenna 
elements. However, they do not include coupling effects due to sub- 
structure contributions which return back to the antenna aperture. 
One way to include these effects is by re-executing the FEM code 
with the secondary external diffraction contribution as part of the 
excitation in addition to the feed. Alternatively, reciprocity can be 
used to first compute the secondary fields arriving at the antenna 
aperture due to plane wave incidence before executing the FEM 
code. 

APATCH Code: This code was developed at DEMACO [24] and 
models the aircraft structure using facets. It is an outgrowth of the 
XPATCH code [25] and since the facets can be as small as 1/10 of 
a wavelength, it permits a good geometrical fidelity of realistic air- 
craft. The field calculations are performed using PO/PTD for the 
specular and first order diffraction contributions. An advantage of 
using the PO is the inclusion of contributions from coated surfaces 
using the plane wave reflection coefficients of the multilayer coating. 
Another advantage of the code is the inclusion of multiple specular 
interactions using the SBR method [26]. Shadowing is typically a 
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major task when hundreds of thousands of facet elements are used 
for modeling the structure, but can be rapidly executed by making 
use of the Z-buffer of the SGI workstation. 

The APATCH code accepts as input either the surface magnetic 
current elements computed by the FEM code over the entire non- 
metallic aperture of the printed antenna or the corresponding volu- 
metric radiation pattern. If the source information is given by the 
magnetic currents, the code generates ray bundles which are then 
interacted with the aircraft structure using the principles of ray and 
physical optics. Each element is treated individually and it is there- 
fore appropriate to group nearby elements to reduce CPU time. The 
final pattern is computed by a linear addition of the contributions 
from all current elements or groups of elements. If the APATCH 
source is the overall antenna pattern in the absence of substruc- 
ture interactions, the code introduces weighted ray fields which are 
subsequently interacted with the aircraft structure using ray and 
physical optics principles. There is enough flexibility into the code 
to compute the overall pattern by adding the primary and secondary 
patterns each provided by different procedures. For example, the fi- 
nal pattern may be computed as the linear addition of the primary 
pattern given by the FEM code and the secondary pattern due to 
the interference of the substructure as provided by APATCH. 

The APATCH code includes a well developed graphical interface for 
rendering the entire aircraft geometry using the open GL library of 
the SGI workstation. Ray interactions can also be displayed along 
with the surface field strengths and this is particularly useful in 
identifying the major disturbances due to the substructure. 


5 Comparison of Measurements and Calculations 

To examine the effectiveness of the proposed hybridization of fi- 
nite element and high frequency codes, the cylinder+wing structure 
shown in Figure 5 was constructed by Naval Air Warfare Center 
Weapons Division (China Lake, CA) 1 and the University of Michi- 

1 R. Sliva and H. Wang provided the original cylinder model with the antenna cavities. 
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gan. The patch antenna geometry is shown in Figure 13 and was 
placed on the cylinder’s surface at a = 28.7°, 45°, and 90° (see 
Figure 5). The radiation patterns of the patch antenna in the pres- 
ence of the wing were obtained at the facilities of Mission Research 
Corporation (Dayton, Ohio) and are shown in Figure 14. It is seen 
that the patterns at a = 28.7° and a = 45° are affected quite 
substantially by the presence of the plate since the plate is visible 
to the antenna. However, when the patch is on the cylinder’s top 
(a = 90°), the plate’s effect diminishes and the overall pattern is 
nearly identical to that in the absence of the plate. 

The computations for the configuration in Figure 5 were done by 
first using the FEMA-CYL code to obtain the surface electric fields 
over the non-metallic portion of the aperture housing the patch. 
These were then grouped in patches of 3x3 pixels in size and turned 
into equivalent magnetic current infinitesimal dipoles whose stength 
was set equal to ( zM z + a 8(f>8z, where M z and are the 

surface current densities over the pixel and a8<f>8z is its area. The 
calculated radiation patterns for the case where a = 45° are shown 
in Figure 15 and seen to be in good agreement with measurements. 
In particular, the pattern based on the UTD code is in better agree- 
ment with the measured data in the shadow region (below the plate) 
but both high frequency codes give nearly identical results in the lit 
region. This is primarily because APATCH employs the less accu- 
rate PO integration to obtain the fields in the shadow region below 
the plate. In contrast, the UTD code uses geometrical optics and 
diffraction theory which is known to be more accurate for flat plates 
with straight edges. This better performance of the UTD code is 
even more pronounced when the patch antenna is brought closer to 
the plate’s surface as shown in Figure 16. Clearly, this is due to the 
proximity of the patch to the plate’s surface, resulting in stronger 
secondary currents on plate’s surface whose PO/PTD approxima- 
tion is much less accurate than that provided by the UTD which 
accounts for ray curvature and surface diffraction effects. The UTD 
pattern is in good agreement with the measured data except near the 
shadow boundary associated with the plate. For practical aircraft 
configurations, the antenna will be placed at far distances from the 
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wing and in that case the APATCH code is more attractive because 
of its geometrical adaptability and capability to handle materials in 
the context of the PO approximation. 


6 Conclusions 

We presented a rather simple hybridization of finite element and 
high frequency codes for antenna pattern calculations in the pres- 
ence of a complex structure such as an aircraft. Basically, the fi- 
nite element code was employed to generate the aperture fields on 
the antennas surface and these were then turned into equivalent 
magnetic currents. These currents were subsequently used as the 
sources to the high frequency codes and the antenna pattern was 
calculated as the sum of the direct antenna radiation pattern and 
that generated by ray interactions with the substructure. The ac- 
curacy of this procedure was verified by comparing the calculated 
results with measured data. Since the measured set-up consisted of 
canonical components, the hybridization with the UTD code leads 
to more accurate results in the shadow region. However, for general 
airframe configurations, the hybridization with the APATCH code 
which combines the PO, PTD and SBR methods is more attractive 
due to its geometrical adaptability and capability to handle material 
coatings. 

Although we did not consider antenna loading effects due to the 
substructure re-radiation, this can be easily incorporated into the 
analysis by executing the finite element code in the presence of the 
antenna feed excitation and the secondary excitation arriving at the 
antenna aperture after undergoing reflection and/or diffraction. 
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Figure 1. Triangular surface mesh of ridge horn used for a volume 
tetrahedral mesh. Note that the ridge is not in the principal plane 
of the horn. 
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Figure 2. Some commonly used tessellation elements for generat- 
ing volume mesh. 
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Figure 3. Printed antenna configuration recessed in a treated metal- 
lic surface. 
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Figure 4. Illustration of antenna pattern illumination and inter- 
action with platform. 
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Figure 5. Measurement model for evaluating the hybrid FEM-High 
Frequency implementation. 



Figure 6. (a) Antenna array on cylindrical platform, (b) FE-ABC 
modeling of an antenna element 
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Figure 7. Comparison of measured and calculated radiation pat- 
terns for a rectangular patch on a cylindrical platform, (a) without 
dielectric covering, (b) with dielectric covering, (c) test body: cir- 
cular ogive cylinder. 
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Figure 8. Input impedance calculations using an improved coax- 
ial cable modeling. 



Figure 9. Modeling of patch antennas using tetrahedral elements. 
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Figure 10. Illustration of an aperture fed dual circular patch an- 
tenna. The portion of the cavity above the slot is modeled using 
tetrahedrals whereas the lower microstrip feed is modeled using brick 
elements. 
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Figure 12. Accuracy of the artificial absorber truncation schemes 
for input impedance computations. 
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Figure 13. Patch antenna geometry and material properties. 
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Figure 14. Measured pattern of the patch for the configuration 
shown in Fig. 5. The shown pattern in the absence of the plate was 
calculated using FEMA-CYL 
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Figure 15. Comparison of measured and calculated patterns when 
the pattch antenna is placed at a = 45° using the hybridization of 
the finite element and high freguency codes. 
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Figure 16. Comparison of measured and calculated patterns when 
the patch antenna is placed at a = 28.7° using the hybridization of 
the finite element and high frequency codes. 
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